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Abstract 

We derive the nonlinear equations satisfied by the coefRcients of linear com- 
binations that maximize their skewness when their variance is constrained to 
take a specific value. In order to numerically solve these nonlinear equations 
we develop a gradient- type flow that preserves the constraint. In combination 
with the Karhunen-Loeve decomposition this leads to a set of orthogonal modes 
with maximal skewness. For illustration purposes we apply these techniques 
to atmospheric data; in this case the maximal-skewness modes correspond to 
strongly localized atmospheric flows. We show how these ideas can be extended, 
for example to maximal-flatness modes. 
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1. Introduction 

When dealing with large data sets it is c onvenient and co stumary to make 



use of the Karhunen-Loeve decomposition (jKarhunenl . 119471 ) in order to ex- 
press the data in terms of the so-called empirical orthogonal functions (EOFs), 
these are linear combinations of the original variables whose second-order cross- 
correlations vanish, i.e. they are linearly uncorrelated, as in ([5]). A typical 
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application consists in reducing then the number of degrees of freedom to a rel- 
atively small number of EOFs with a variance larger than a certain threshold, 



for more applications see for example (jPreisendorfer 



19881 ). In many systems 



of interest the EOFs are nonlinearly correlated, a fact that can have important 
consequences. Similarly, the probability densities of the time-dependent EOFs' 
amplitudes are often approximately Gaussian but the deviations from Gaus- 
sianity may be of great relevance. One indicator of nonlinear correlations and 
of the non-Gaussian character of fluctuations is the skewness, the third-order 
moments of the variables' distribution. In this article we show how to construct 
orthogonal linear combinations of the variables that maximize the skewness and 
we present a numerical method in order to solve the ensuing nonlinear equa- 
tions. For the purpose of illustration we apply it to meteorological data; the 
maximal-skewness modes found in this case correspond to spatially localized 
and meteorologically meaningfuU atmospheric flows. 

2. Maximal skewness modes 

Given a set of n dynamical variables, 

Vi{t), 1 < i < n, 
with vanishing mean we want to construct a linear combination iplt) 

n 
i=l 

such that its skewness s, 

s(ai, . . . ,a„) := (■0^) , 

is maximal. The angular brackets indicate an average over the data set. In 
order to get a sensical solution some restriction must be imposed on it. An 
appropiate restriction is to fix its variance, say 

(V^^) = 1. (1) 
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This choice imphes that the dimension of is [vi] . In terms of the coefficients 
{tti} the constraint reads, 

n 

^ d^C^^ci, ^ 1, (2) 

where are the elements of the covariance matrix Cu := (viVi) . The maximal 
skewness is obtained by setting equal to zero the variation of (J — A^V^)) 
where A is the Langrange multiplier associated with the variance constraint (fT|). 
In this way one obtains 

n n 

3 ^ afcS'feHfl; - 2A^ Cija; = 0, I < i < n, (3) 

k,l=l 1=1 

with 

Skii := {vkViVi) , 

the elements of the skewness tensor S. We assume all the elements of the tensors 
C and 5 to exist. The values of the n coefhcients Ui and of the Lagrange 
multiplier A satisfying the (n + 1) equations (U) and ([3]) will be denoted as 
af and as A^ respectively. The number of real solutions maybe larger than 
n, see for example Fig. [TJ If {Xa,af} is a solution then also {— A^,— af} is 
a solution. The value A^ associated with a solution a" is proportional to the 
skewness s^(a", . . . , a"). This can be seen by multiplying the i-th equation ^ 
by ai summing over all the components and making use of the constraint ([2|, 
one obtains then that 

2A«=3s(a?,...,af). (4) 

Without loss of generality, we can take the {vi{t)} to be EOFs, i.e. their 
covariance matrix Ca := (viVi) is 

Cii = WiSa- (5) 
Then the equations to solve become, 

3 ^ akSkuai ~ 2\wfa.i = 0, I <i <n. 

k,l 
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It is convenient to introduce the dimensionless quantities 



and S, 



Pi := WiQi, 

{VkVlVi) 



kli 



WkWlWi 

In terms of these equation ^ reads, 



3 ^ PkSkuPi - 2Aft = 0, 1 < z < n, (6) 

k,l 

and the unity covariance constraint ^ is, 

n 
i 

A more compact way of writing equation ^ is 

^(p) = 2A?, 

where It is the gradient of the skewness function s(/3i, . . . , :— ?(ai, . . . , a„), 
i.e. 

^.(^):=^ = 3E5.,fe/3,/3fc. (8) 

The solutions to equations ( (5]) and ([7]) may correspond to saddle points of 
s(/?i, . . . this is further analyzed at the end of Section [3l With n ~ 2 the 
solutions can be expressed analytically, see the Appendix. With n > 2 one has 
to find them numerically, in the next Section we present a way of doing this. 

If all the Wi 's have the same dimensionality then we can associate a variance 
to each solution namely 

n 

W^:=Y,m'wl (9) 

i 

It may happen that a solution has a relatively large skewness s{/3f, . . . , (3") while 
its variance is relatively small. Usually one is interested in solutions with 
both quantities relatively large. Accordingly one can consider a dimensional 
skewness parameter, call it 

E,.:=s(/3i",...,/3^)VF3. (10) 
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In closing, let us mention that we take the space of the coefficients {(3i ,...,/?„} 
to be Euclidean, i.e. the inner product of two /3— vectors is 

y n 

and the constraint "^Pf = ^ means that the vectors j3 are of length 1. 

3. An algorithm in order to solve the system of equations ([6]) and ([7]) 

Consider a gradient flow, i.e. let (3 (r) move downhill a potential V{ (3)^ 

dp^ _ dV 
dr d(3i ' 



so that V{[3) is a Lyapunov function, 

. dV d(3. 



dy_ 

dr 



E 



dPi dr 



dV 

dl3. 



< 0, 



and the evolution of f3 (r) stops when an extremum of V{P) is reached. In 

» 2 



In this context notice 



general, such an evolution will not conserve 
that if {Ac/Jf} is a solution of ([6]) then also {/iAa, /i/?,f } is a solution of ([6|). 
Therefore, we can work with /^-vectors of arbitrary length if we replace A by 



A, i.e. instead of (El) we can solve 



^(/3) = 2A 



(3. 



f3 



(11) 



in this formulation the 



Since both the Ihs and the rhs are proportional to 
length of P does not play any role. Once a vector satisfying this equation is 
found, then we normalize its length in order to obtain a solution of the equations 
® and 0. 

The previous considerations lead us to introduce the following deviation 
vector A ( /3 ) , 



A(/?):= 



"?(/?) -2A 

— > 2 



P 



P 
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and a potential V{P) given by 
» 2 



(12) 



12A 



/3 



s(/3) + 4A2 > 0. 



The deviation vector A ( /3 ) does not depend upon the length of /3 , it vanishes 
when P solves (jlll) and A equals the corresponding Lagrange multiplier Aq. 
Therefore, the value ofV{f3) measures the departure of {A, /? } from a solution 
to (jlip and, after normalization, from a solution to equations (|6|7p . This specific 
form of the potential has been chosen because, as it will be seen, it conserves 



the length 



. From it one finds that, 



dT 



dV 



41^1 



36As 



/3 



f3 



^3^ 



12A 12 ^ 

/3 kj 



f3 



(13) 



One can check that indeed this implies d 



/dr = 0. Therefore, we can fix 



/3(r) 



= 1 and get 
dr 



with 



l/3| = l 



[4a2 - 36As] A + 12Aa, - 12 ^ S,,kPjak (14) 

kj 



=0. 



(15) 



In the definition of the potential V{f3), equation (fT2)) . and in all the equa- 
tions we have derived from it, A appears as a free parameter. In order to find 
the solutions to the equations ([6]) and ([7]) A must take the value Aq, associated 
with the vector solution /?". Therefore one should let also the Lagrange multi- 
plier A evolve until it reaches the searched value Aq. This is achieved by taking 
A(r) — \m{P) where \m{ P ) is the A value that minimizes V{(i), 

= 0^8Am Z?" 

A— Am 



dV_ 



1.6. 



s{t). 



12s = 0, 
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Compare this with At this value of A the potential V{/3) equals, 



Vm{P) V{(3) 
-i 



A=Aj 



^{13) 



(16) 



The evolution of the (3 (r) induced by this potential is 



dT 



+ - 
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/3 fcj 



(17) 



In agreement with the results in the previous paragraph, also these equations 



lead to d 



f3 



/dr = 0. Therefore, one could simplify them by taking /3 (t) 
1, however, in order to control numerical errors, it is preferable to use these 
equations as they stand. As before, one has that 



dVi 



M 



dT 



E 



dVM dl3, 
dj3i dT 



-E 



M 



dp. 



< 0, 



i.e. Vm{ /? ) is the corresponding Lyapunov function. 

In Fig. [T] we see an example with n — S constructed from the meteorolog- 
ical data used in Section Fig. [1] shows the level plots of s{(32, Pt, Pg) and 
of the corresponding potential V((32, Pt, Pg) on one hemisphere of the sphere 
/3| + Z?! + = 1. The angle ai = arccos(/32) and the angle a2 = arctan(/37//3g). 
The blacks dots indicate the positions where the potential vanishes, these co- 
incide with the positions of the two maxima, two minima and three saddles 
of s(/?2, /?7, /Jg). In addition. Fig. [T] shows how, starting from initial values of 
P2, Pt and Pg, on a regular lattice the gradient-flow equations ^T7\ lead to these 
critical values. Whether a solution /3" is an extremum or a saddle depends upon 
the character of [hfA the Hessian matrix at the solution projected on the sphere 



P 



1 , i.e. by the character of the {n — 1) x (n — 1) matrix with elements 



6°' 



ip^r 



with 1 < i,j <n — 1, 



7 



where 

n 

:=6 5^/3^5fcH l<l,i<n, 

k 

is the nx n Hessian at the solution /3" and cr" :— (Tj (/?"). In these expressions 
it is assumed that ^0. 

4. Numerical implementation of the gradient algorithm 

In the general case one has that the number of extrema grows explosively 
with n. We deal with this problem by taking first the ten largest EOFs and 
using pT]) with 1000 random initial values of . . . ,(3io} in order to find the 
combination with the largest skewness, then we increase the number of EOFs 
by taking the fifteen largest EOFs and as initial /3-values the solution found in 
the previous step suplemented by /3ii ==••• = f3i^ — and let the /3's evolve 
according to (jl7p with noise added to them. Once a new solution is found the 
last step is repeated but now with twenty EOFs, etc. The noise is added in 
order to explore larger sections of the phase space and not to remain trapped 
in a local extremum. A 4th-order Runge-Kutta algorithm is used in order to 
integrate the system of ODEs (|17l) . A solution is found when the value of the 
potential V is close enough to zero and the value of the skewness has virtually 
ceased to change. 

In Fig. [2] we show the results obtained when this procedure is applied to the 
meteorological data used in Section ([5]). One can see an increasing maximal- 
skewness value as n increases from 10 to 50. 

5. Orthogonal set of maximal-skewness modes 

In order to generate a set of linearly uncorrelated orthogonal modes ordered 
according to their skewness one should proceed as follows. Firstly, the method 
presented in the previous sections is used in order to create from the data 
m{x,t), 1 < a; < n, the linear combination with the largest skewness, which we 
write as follows, 

n n 
i=l i=l 
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i.e. the ei{t) are the normalized EOFs, 

ei{t) ^ w:[^Vi{t), with (viVi) ^ w'^dii \<i,l<n. 



Thanks to the Karhunen-Loeve theorem ( Karhunen 



1947 ) each Vi is associated 



with a unique n-dimensional eigenvector 7ri(x) satisfying 

n 
x = l 

= {m{y, t)m{x, t)) , l<x,y<n. 
Therefore, if 7^ w^^, the corresponding eigenvectors are orthogonal, 

71 

x=l 

and can be taken to be normaUzed. If the discrete indices x and y are associated 
with positions in physical space then each of the eigenvectors TTi{x) describes 
a spatial pattern. In such a case, to ip^{t) = X]"=i '^i^i(^) there corresponds a 
unique spatial pattern 'iIj^{x) :— X]r=i 

Using the covariance metric ([5|) this V'^(i) rnode is projected out from the 
n-dimensional set {ei(t), e2(t), . . . , e„(t)}. The new set so obtained is 

m,{t) = (1 - f3l)e,{t) - A ^ Pueu (t) , 

{'ip^{t)m,{t)) =0, 1 <i<n. 
Their covariance matrix is 

{rhj{t)mi{t)) = -PiPj, i ^ j. 

It has one vanishing eigenvalue with eigenvector ip^it) = X]"=i Pi^ii^) (n — 
l)-times degenerate eigenvalue 1 with normalized eigenvectors 

Skit) := {^l+f3iy^ i-pkmiit) + Pimkit)] 
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By construction the efc(i)'s and ip^{t) are linearly uncorrclated (^^jj^{t)ek{t)') = 0. 
To each ek{t) there corresponds a unique pattern efc(x) oc — afc7ri(x) + aiT:k{x). 
It follows that ?/'^(a;) and the (n — 1) patterns ek{x) are orthogonal since 

n 

^^il^^{x)ek{x) oc -afcfli +aiafc = 0. 

x=l 

As indicated in (fTUl) instead of V'^(i) one may consider the associated dimen- 
sional mode Wiip^{t). 

6. An application to meteorological data 

Using the data available at the ECWMF ERA40 website we computed 
the daily values of the streamfunction anomalies on the Northern-hemisphere 
500 hPa level during the months December, January and Febr uary from 195 8 



to 2001. The skewness in t his field was already noticed by 



White! dlQSOl) : 



Nakamura and Wallace! ( 19911 ). The dimensionless skewness field s (a:) := {rn?{x^t)) (m^(a;,i)) 
is shown in Fig. [3l Two black dots indicate the two positions with the largest 
positive and negative skewness s(a;), to wit: (168 East, 47 North) with skewness 
0.76 and (158 East, 19 North) with skewness -0.39. 

Fig. d] shows the spatial patterns corresponding to ^jj^ the largest and ijP' 
the second-largest skewness mode. The skewness of these modes are si — 1.13 
and S2 — 0.96. That these values are larger than the maximal values of the 
local skewness s{x) indicated by the black dots in Fig. [3] is possible due to the 
non-vanishing third-order cross-correlations {m{y,t)m{x,t)m(z,t)) ^ for x,y 
and z not simultaneously identical. 

For these data one has that X]i=i ~ ^■^J2i=i '^h i-^- the first 50 EOFs 
describe 90% of the total variance of the daily streamfunction fields. In par- 
ticular, wf and w| describe 7.7% and 6.4% of the total variance respectively. 
One also finds that the variances Wi and corresponding to the maximal- 
skewness modes ^/i^ and V'^i confer account for 3.3% and 3.9% of the total 
variance. Their contributions to the total variance is comparable to those of 



-3/2 
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the EOFs vg and respectively. Therefore, these maximal-skewness modes are 
physically relevant. 



7. Possible generalizations 



For example, instead of the modes with maximal skewness one could be 
interested in the modes with maximal flatness /(ai, . . . , a„) :— (V"^) given that 
(ip"^) = 1. The equation analogous to ^ is 



7] 



k,l.m—l 



2XP, = 0, l<i<n, 

{VkVlVmVi) 



WkWiWmWi 

The Lagrange multiplier A is proportional now to the flatness of the solutions, 
= 2/(a", . . . , a"). The corresponding deviation vector and potential are 



A4(/3):= 



0(/3)-2A 



/3 



/3 



A. 



with 



mi ■ 

k,l,m—l 

The minimum value of ) is achieved when A equals Xm — 2 
the potential is then 



/9 



-6 



0(/3) 



16 



/3 



/'(/3i,...,/3«). 



As one can check J2 Pi {dV^M /dfii) — 0, so that dp /dr = 0. 
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Appendix 

With two EOFs and after eliminating A there is only one equation to solve, 
namely 

/?2 (PlSui + 013211 + 2/3i/325i2i) = /3i (/3?5n2 + /3|5222 + 2/3i/J25i22) • 

Dividing both sides of this cquahty by 01 one gets that Z2 '■= 02 /Pi is the 
solution of the following third-order polynomial, 

Z2 {Sui + Z2S22I + 2^2S'i2l) = (<S'll2 + 2^2'^222 + 222'S'l22) • 

Once Z2 is known, we recover /3i and /32 from 

P2 +f3l = l^ pf^l + = 1 ^ = (1 + ^2^-1 

and = zlil + zl)-\ 
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Figure 1: The black contours correspond to isolevels of the skewness s{f32, /Jy, /3g). The shading 
denotes the values of the potential V {(32 , 07 , l3g) . The black dots indicate the location of the 
skewness maxima, minima and saddles; the non-negative potential vanishes on these locations. 
The white lines are the trajectories generated by equation II17I I with initial conditions on a 
regular lattice. 
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Figure 2: Skewness of the maximal-skcwness mode as a function of the number of EOFs 
included in the calculations. 
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Figure 3: Skewness s{x) of the daily streamfunction on the 500 hPa surface in the months 
December through February, years 1958 through 2002. The black points indicate the positions 
with maximal positive skewness s{x) = 0.79 and maximal negative skewness s{x) = —0.39. 
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Figure 4: The spatial patterns of two maximal-skewness modes of the daily streamfunction 
on the 500 hPa surface in the months December through February, years 1958 through 2002. 
On the left, the largest-skewness mode with skewness si = 1.13 and, on the right, the mode 
with next-largest skewness S2 = 0.96. 
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